int gen_hilbert_recur(int *data, int dim, int n, int s,
		int x, int y, int z,
		int dx0, int dy0, int dz0,
		int dx1, int dy1, int dz1,
		int dx2, int dy2, int dz2){
  if (s == 1) {
    data[(x * dim + y) * dim + z] = n;
    return n + 1;
  } else {
    int sl = s >> 1;
    if (dx0<0) x-=sl*dx0;
    if (dy0<0) y-=sl*dy0;
    if (dz0<0) z-=sl*dz0;
    if (dx1<0) x-=sl*dx1;
    if (dy1<0) y-=sl*dy1;
    if (dz1<0) z-=sl*dz1;
    if (dx2<0) x-=sl*dx2;
    if (dy2<0) y-=sl*dy2;
    if (dz2<0) z-=sl*dz2;
    //int dx0, dy0, dz0, dx1, dy1, dz1, dx2, dy2, dz2;
    int nnew = n;
    //nnew = gen_hilbert(data, dim, nnew, sl,
    // int sl;
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x, y, z, dx1, dy1, dz1, dx2, dy2, dz2, dx0, dy0, dz0);
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x+sl*dx0, y+sl*dy0, z+sl*dz0, dx2, dy2, dz2, dx0, dy0, dz0, dx1, dy1, dz1);
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x+sl*dx0+sl*dx1, y+sl*dy0+sl*dy1, z+sl*dz0+sl*dz1, dx2, dy2, dz2, dx0, dy0, dz0, dx1, dy1, dz1);
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x+sl*dx1, y+sl*dy1, z+sl*dz1, -dx0, -dy0, -dz0, -dx1, -dy1, -dz1, dx2, dy2, dz2);
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x+sl*dx1+sl*dx2, y+sl*dy1+sl*dy2, z+sl*dz1+sl*dz2, -dx0, -dy0, -dz0, -dx1, -dy1, -dz1, dx2, dy2, dz2);
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x+sl*dx0+sl*dx1+sl*dx2, y+sl*dy0+sl*dy1+sl*dy2, z+sl*dz0+sl*dz1+sl*dz2, -dx2, -dy2, -dz2, dx0, dy0, dz0, -dx1, -dy1, -dz1);
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x+sl*dx0+sl*dx2, y+sl*dy0+sl*dy2, z+sl*dz0+sl*dz2, -dx2, -dy2, -dz2, dx0, dy0, dz0, -dx1, -dy1, -dz1);
    nnew = gen_hilbert_recur(data, dim, nnew, sl, x+sl*dx2, y+sl*dy2, z+sl*dz2, dx1, dy1, dz1, -dx2, -dy2, -dz2, -dx0, -dy0, -dz0);
    return nnew;
  }
}
void gen_hilbert(int *data, int dim){
  gen_hilbert_recur(data, dim, 0, dim, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1);
}